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Abstract. In this paper, we propose a solution for a fundamental problem 
in computational harmonic analysis, namely, the construction of a multireso- 
lution analysis with directional components. We will do so by constructing 
subdivision schemes which provide a means to incorporate directionality into 
the data and thus the limit function. We develop a new type of non-stationary 
bivariate subdivision schemes, which allow to adapt the subdivision process 
depending on directionality constraints during its performance, and we derive 
a complete characterization of those masks for which these adaptive direc- 
tional subdivision schemes converge. In addition, we present several numerical 
examples to illustrate how this scheme works. Secondly, we describe a fast 
decomposition associated with a sparse directional representation system for 
two dimensional data, where we focus on the recently introduced sparse direc- 
tional representation system of shearlets. In fact, we show that the introduced 
adaptive directional subdivision schemes can be used as a framework for deriv- 
ing a shearlet multiresolution analysis with finitely supported filters, thereby 
leading to a fast shearlet decomposition. 



1. Introduction 

Efficient and economical representations of anisotropic structures are essential 
in various areas in applied mathematics. The nature of the problems we face can 
be divided into two types, namely when the anisotropic structure is given explicitly 
and when it is given implicitly. The analysis of images and higher dimensional 
data with respect to directional features shall serve as an example of an explicitly 
given anisotropic structure, whereas the solution of hyperbolic partial differential 
equations often exhibits the phenomenon of shocks which can be interpreted as an 
implicit anisotropic structure. 

It is well known that wavelets are perfectly suited for providing efficient represen- 
tations in the sense of sparsity for problems with a dominant isotropic regularity, 
at the same time being associated with a multiresolution analysis which is the 
key ingredient for a fast decomposition algorithm. However, when dealing with 
anisotropic phenomena wavelets do not perform equally well. In fact, it can be 
proven that wavelets do not provide optimally sparse representations. 

In contrast to earlier approaches such as directional wavelets complex wave- 
lets [21], ridgelets [3], and contourlets 16J, the curvelets introduced by Candes and 
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Donoho precisely satisfy this need, in the sense of resolving the wavefront set |5] and 
the curvelet representation being optimally sparse for objects with C^-singularities 
Also there already exist some first results on applying curvelets to hyperbolic 
partial differential equations by Candes and Demanet [2 . However, one drawback 
is the lack of a multiresolution analysis associated with curvelets, and, in particu- 
lar, a fast decomposition algorithm in the time domain. This raises the question 
about the existence of a representation system with analyzing properties as good 
as curvelets, but being equipped with a more "wavelet-like" structure in the sense 
of being associated with a multiresolution analysis. In fact, the discrete counter- 
part would then lead to finitely supported filters that allow for a mathematically 
justified discrete fast decomposition of discrete data. We anticipate such a represen- 
tation to combine the favorable computational properties of wavelets with the main 
additional property to provide a means to resolve anisotropic structures efficiently. 

In this paper we give a complete, positive answer to the question of the existence 
of such a system by introducing subdivision schemes for the recently introduced con- 
cept of shearlets, thus constructing an associated multiresolution analysis which 
indeed leads to a fast discrete decomposition algorithm. The directional represen- 
tation system of shearlets [IQJ stands out for the following reason. They do not only 
precisely resolve the wavefront set [25] and provide optimally sparse representations 
PU] , but shearlet systems are generated by one single function which is dilated by a 
parabolic scaling and a shear matrix and translated in the time domain, hence form 
an affine system. We might even interpret the system of shearlets as being gen- 
erated by a strongly continuous, irreducible, square-integrable representation of a 
certain group, the shearlet group TT\. This rich mathematical structure enables, for 
instance, the application of coorbit theory to study smoothness spaces - so-called 
shearlet coorbit spaces - associated with the decay of the shearlet coefficients [T^ . 
We would further like to mention that one attempt to associate shearlets with a 
so-called generalized multiresolution analysis can be found in [24 . However, this 
structure did not yield a fast decomposition due to the fact that the filters are not 
compactly supported and even infinitely many filters have to be employed. 

Our approach to derive a multiresolution analysis associated with shearlets and 
to provide a feasible fast shearlet decomposition comprises the introduction of a 
new class of non-stationary bivariate subdivision schemes which incorporate direc- 
tionality in a particular way. Subdivision schemes provide a mathematical method 
to refine given coarse data while providing characterization results to ensure conver- 
gence to a continuous function, say. Moreover, such schemes automatically provide 
refinable functions which are the basis for any multiresolution analysis as nest- 
edness of the different levels of resolution is equivalent to the refinability of the 
underlying "basis" function. Homogeneous stationary subdivision schemes have 
been studied extensively over the last 20 years; for an elaborate survey we refer 
the reader to [5]. Recently, algebraic methods have been introduced as a means to 
derive characterizations of convergence and approximation order in a very natural 
way for multivariate subdivision (cf. [29]). On the other hand, also the conditions 
of homogeneity and stationarity have been released by various authors, leading to 
subdivision schemes where the refinement rule varies with the level of iteration or 
the location of refinement. However, the gain in generality always comes with the 
prize of a loss of structure so that there is comparatively little known about these 
generalizations (see, e.g., [9l|7]). In particular, no subdivision schemes were known 
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SO far which provide a means to adapt the subdivision process depending on direc- 
tionality constraints during its performance while still ensuring convergence. The 
development of such subdivision schemes will be important both for construction 
of a shearlet multiresolution analysis as well as for opening the research area of 
methods for data refinement to incorporate anisotropic structures. 

We will show in this paper that such an adaptive directional subdivision scheme 
can be constructed and it will indeed lead to a shearlet multiresolution analysis 
and a fast shearlet decomposition. Our approach to derive a non-stationary bi- 
variate adaptive directional subdivision scheme is based on the idea to iteratively 
apply two subdivision schemes each of which is associated with a different direction. 
The two individual subdivision schemes can employ two different finitely supported 
filters while their respective dilation matrices are taken from the theory of shear- 
let systems. We would also like to mention at this point that the most natural 
"directional" operation, the rotation, can not be employed, since its action does 
not provide a refinement of a lattice. In contract to this observation, products 
of parabolic scaling and shear matrices do indeed satisfy this desirable property. 
The constructed subdivision scheme provides the opportunity to adaptively change 
the orientation of the data during the subdivision process, since in each iteration 
one of both single subdivision schemes can be applied. In this sense, we can visu- 
alize the subdivision process as a binary tree, in which the direction of the finer 
data is dependent on the branch we choose. However, for convergence we certainly 
need to study each branch of the tree, which requires an appropriate definition of 
convergence. Our first key result shows that, provided the adaptive directional sub- 
division scheme converges, we obtain associated generalized refinement equations 



(Theorem 4.6 1. These will become essential for deriving a shearlet multiresolution 
analysis. As a main result we then provide a complete characterization of those 
masks which lead to convergent adaptive directional subdivision schemes (Theorem 



4.141 in terms of algebraic and spectral properties of the associated filters. In the 
proof we will make use of ideal theoretic methods which come in handy to extract 
"the zero at —1" of the two masks. 

For the construction of a shearlet multiresolution analysis we employ the fact 
that each wavelet multiresolution analysis is associated with a convergent subdivi- 
sion scheme [T4]. We introduce scaling spaces based on the previously constructed 
directional subdivision schemes, and then prove that these indeed provide a mul- 



tiresolution analysis structure (Theorem 6.3) due to the refinement equations men 



tioned above. This multiresolution analysis will then provide us in a very natural 
way with a mathematically justified discrete fast shearlet decomposition of discrete 
data which is stated as Algorithm |7.6| Also here we encounter a binary tree struc- 
ture, since the decomposition will be dependent on the different directions which 
were encoded in a binary tree structure of the subdivision process. For the con- 
struction of a shearlet multiresolution analysis and a fast shearlet decomposition, 
we focus on the situation of interpolatory masks. The non-interpolatory case is 
beyond the scope of this paper and will be studied in a forthcoming paper. 

The outline of the paper is the following. In Section [2] we briefly introduce 
discrete shearlet systems. We further study which directions can be attained by 
the action of the associated dilation matrices on 1? . The new type of subdivision 
schemes, which we baptize adaptive directional subdivision schemes, are introduced 
in Section [3] In Section |4] we provide a complete characterization of convergence for 
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those schemes along with the necessary ideal theoretic background. Some numeri- 
cal experiments on the refinement of data employing this new type of subdivision 
schemes are provided in Section |5] We then show how the previously derived 
adaptive directional subdivision schemes can be used as a framework for deriving 
a shearlet multiresolution analysis with finitely supported filters (Section [6|. In 
Section [7] we employ these results to provide a fast shearlet decomposition. 

2. Refinement of by Anisotropic Scaling and Shearing 

2.1. Shearlet Dilation Matrices. Our approach towards directional refinement 
of the lattice and, later on, adaptive directional subdivision schemes is inspired 
by the recently introduced discrete shearlet transform |19) . since this transform 
is able to precisely detect directions of singularities (cf. [12]) which we will take 
advantage of. In order to provide a thorough motivation for our construction, allow 
us to first briefiy review the idea of shearlets. 

Each shearlet system forms an affine system, i.e., consists of dilations and transla- 
tions of one single generating function t/j £ _L^(]R^), a so-called shearlet. As dilation 
matrices, products of anisotropic parabolic scaling matrices and shear matrices - 
which coined the name "shearlets" ~ are employed. In order to define a shearlet 
system, let Aa, a > 0, and Ss, s G M, which are defined by 

^a-(o ^) and Ss^l^Q 7), 

denote a parabolic scaling matrix and a shear matrix, respectively. Then the shearlet 
system associated with a shearlet tp e L^(IR^) is given by 

(2.1) {i^jkmix) 2-iX5_fc^4-.2; - m) : 3, k e Z, m £ 1?}. 

The three parameters j,k,m are interpreted in the following way: j provides the 
scale, and k and m detect the direction and position of singularities, respectively. It 
is easy to construct shearlets such that ( |2.1[ ) forms a Parseval frame for _L^(IR^), for 
instance, by choosing ^'(^1,^2) = fpi{£,i)i^2{£,2/S,i) , where V'l G L'^{R) is a discrete 
wavelet, i.e., X^^gz IV'i(4"''^)P = 1 for w € M, satisfying ^pi £ C°°(M) and supp?/;i C 
[-1,-i] U [i,l], and e L'^{R) is a bump function satisfying ■02 € C°°(M), 
supp ^^2 C [—1, 1], and J2kez l^^ik + uj)\'^ — 1 for uj £ R (cf. [in])- The associated 
Shearlet Transform 57i^, is then defined on L^(M^) by 

ST-C^fU, k, m) = (/, ijjkm) ■ 

In order to provide an equal treatment of the direction of the x- and y-axis, the 
frequency plane is split into the cone 

C = {{i^,a)£R'■. |ei|>i, |||<1}, 

its by 90° rotated copy, and the square centered at the origin of side length ^ . The 
Shearlet Transform acts on C and its copy as described above, while the choice of 
tp has to be adapted appropriately. The center square can be filled in such a way 
that this system also forms a Parseval frame. The shearlet system in C and its 
copy is usually referred to as shearlets on the cone, see The associated tiling 
of the frequency plane is illustrated in Figure [T] 
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The refinement matrices interesting to us for deriving a directional refinement 
of the lattice 1? are the dilation matrices used in (2.1 1 for j = 1, i.e., the matrices 

Following the philosophy of the shearlets on the cone, also the matrices 

which serve as dilation matrices for the rotated copy of C, will be employed as 
refinement matrices. The matrices and Mk not only provide the possibility to 
map a line to various directions, but moreover possess the property of refining the 
lattice 1? equally at each level as it is shown in the following result. 

Proposition 2.1. The following conditions hold. 

(i) For all j, fc G Z, we have 

Mk e GLaW and Mfc(4^^Z x = A-^^+^^Z x 2-'^^+^^Z. 

(ii) For all j, fc € Z, we have 

Mk G GLaW and Mfc(2"^Z x A'^Z) = 2~^^+^'>Z x 4~(j+'^^Z. 

Proof, (i) The first claim is obvious. To prove the second claim, let j, /c G Z and 
m — {mi, 7712) G Z^. Then 

^2-Jm2y ^ 2-(J+i)m2 ) ^\ 2-0+i)m2 

which implies Mfe(4-JZ x C 4-0+i)z x 2-^J+^^Z. 

Now let n = (711,712) G Z^. Then choosing m = (7711,7712) G Z^ as 7711 
77i — 2^''"^fc77,2 and 777,2 = "•2 yields 

^A-^mA _ /4-(J+i)(7ii - 2J+ifc772 + 2J+iA:772)\ _ A-^+i)^^^ 

,2-^7772^^1^ 2-0 + 1)772 j " 1^2^(^ + 1)772, 

Thus A4(4-JZ X 2-JZ) D 4-(-''+i)Z x 2-0+i)Z, which proves the claim. 



Mk 



Mk 
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(ii) This follows by using similar arguments as in part (i). □ 

Thus, when applying a sequence of matrices M}^^ , ■ • • , iteratively to the 
lattice Z^, at the jth level the points {(4--'(toi +l\),2-^{m2 + \)) : i € {1,2,3}} 
are added to the lattice 4~^Z x 2^-'Z. This is true for an arbitrary choice of integers 
kj g Z, 1 < j < n. Moreover, at each level this map is bijective. 

A similar result holds for the matrices Mj., fc € Z. 

2.2. Feasible Directions. Let us now delve deeper into the explicit construction 
of the refinement by using the splitting idea of the shearlets on the cone. The 
overall aim is to provide a way of refinement such that the points on the y-axis - 
or any other line through the origin - can be moved to an arbitrary line through 
the origin during the refinement process. This immediately forces the refinement 
scheme to provide different strategies for refinement. We will see how this is can 
be achieved by using the matrices and even only for e = — 1,0, 1. In the 
sequel we will only focus on the matrices M^, e = —1, 0, 1, since the others can be 
treated simultaneously. 

In the very first step of the refinement, we apply to 1? for e = —1,0,1. 
Application of e = does not change any directions, £ = 1 maps the y-axis to the 
angle bisector in the first and third quadrant of the plane, and e = —1 has the 
same effect on the second and fourth quadrant. From now on, we consider the two 
cases e G {0,-1} or e € {0,1} separately. Focusing on the second case, in each 
step we not only derive the refinement from a coarser scale A~^'L x 2~^Z to a finer 
scale 4'(J+i)Z x 2-(-'+i)Z, but also have two different ways to achieve this, either 
by applying Mq or by applying Mi. Hence, at the nth level we have applied a 
product of the form M^^ . . . M^^ to Z^, where £j G {0, 1} for each 1 < j < n. For 
e G {0, —1}, one can proceed in exactly the same way which we will, however, not 
work out in detail in this paper. 

From now on, we will use the abbreviation En — {0, 1}", n G N, for the index 
sets and will also denote by 

riGN 

the set of all finite 0-1-sequences and by E^o — {0,1}^ the space of all infinite 
sequences. Note that E is canonically embedded in E^ by the mapping 

E3 e* ^{e,Q,Q,...) &E^. 

The main question to ask at this point concerns the possible directions this 
procedure allows us to map the points on the y-axis to. For this analysis, we 
restrict our attention to the first quadrant of the plane, since the same refinements 
occur in the third quadrant only in an origin-symmetric way. 

We first notice that the sequence of n matrices we choose is completely 
determined by the associated sequence e G En. Hence this refinement scheme has 
the structure of a binary tree as illustrated in Figure [2] 

The directions which might be obtained employing this refinement scheme are 
encoded in this binary tree in a special though natural way. To explore this rela- 
tion, we first compute the product of the matrices which is applied to achieve the 
refinement at level n. Interestingly, the following binary number appears therein. 
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MoMqI?, £=(0,0) 




MqM^I?, £=(1,0) 
MiMoZ2, £ = (0, 1) 



MiMiZ^, £=(1,1) 



Figure 2. The binary tree up to level 2 associated with the re- 
finement scheme. 



Notation 2.2. For e e En, n e N, we define 

(£)2 = X! 2^ Me = • . . . • . 

Using this notion we obtain the following form for a refinement matrix 
Lemma 2.3. Let n e N and e G En. Then we have 

2-" 



Proof. We will prove this lemma by induction. For n = 1, the claim obviously holds. 
Now suppose that the claim is true for some n e N. Let e = (£',£„+i) e En+i, 
e' e En. We have to distinguish between £„+i = 0, hence {e)^ = (£')2' '^^^h 

^4-1 \ A-" 4-"2 (£')2A /4"^"+^^ 5 4-"(£')2' 

2-7 2-" y 1^ 2-("+i) 

Tl+l 



and £„+i = 1, i.e., {e)2 = (£')2 + 2""^ , where 



4-1 2-i\ A-" 4-"2(£')2\ _ /4-("+i) i (4-"(£')2 + 2-") 



2-^ J \ 2-"- J \ 2-("+i) 

4-(n+i) 1 4-n ((e')2 + 2")\ _ /4-("+i) i 4-" (0)2' 
2-("+i) j ~ I 2-("+i) 



which advances the induction hypothesis. 



□ 



Notation 2.4. Let _L be a line through the origin and e G E^, n E N. Then s(L, e) 
denotes the slope of M^L, which is again a line through the origin. We further 
write s{L) for the slope of L. 

The next result computes the values of the slopes s{L,e). 
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Lemma 2.5. Let L be a line through the origin and e G En, n e N. Then the 
following relations between s{L,e), e and the original L hold. 

(i) If L is a line through the origin with s{L) G (0,cxd), then 

(ii) If L ^ {0} X E, i.e., s{L) = oo, then 

2n — 1 

s{L,e) 



(eh 

where we set 2"^^/0 :~ oo. 
(iii) IfL = Rx {0}, I.e., s{L) = 0, then 

s{L,e) = 0. 



Proof, (i) We consider the point (l,s(L)) G L. Using Lemma 2.3 we compute 



1 

s{L) 



r" 4-" 2 {e)2 
2-" 



1 

s{L) 



A--{l + 2s{L) {e)2, 



Hence, the slope of the hne M^L equals 

4"s(L) _ 2"s(L) 

2»(l + 2s(L) (£)2) ~ l + 2s(L) (£)2 " 

(ii) Here we consider the point (0, 1) G L = {0} x 
12.31 we obtain 



s{L) 



2(e)2 

Again employing Lemma 



4-« 4-"2(e)2 
2"" 



4-"2(£)2 

2-n 



Thus 



s{L,e) 



4" 



2" 



2»+i(£)2 2(e)2' 
(iii) is easily verified by noting that the point (1,0) is mapped to 



r" 4-" 2 {e)2 
2-" 



4-" 




so that the slope remains zero. 



□ 



Our main result in this section will show that indeed the points on an arbitrary 
line through the origin of slope 7^ can be moved arbitrarily close to prescribed 
lines through the origin during the refinement process. 

Theorem 2.6. Let L be a line through the origin with s{L) G (0,oo]. Then, for 
each t G 00] and S > 0, there exists some n G N and e G En such that 

\s{L,e)-t\ < S. 

Proof. Suppose L is a line through the origin with s(L) G (0,oo). The case s(L) — 
00 can be dealt with in a similar way. 

For given t G (|, 00) and S > 0, due to the denseness of rational numbers there 
exists some n G N and e G En such that 



n-l 

3=0 



■n+l 



1 



^t{t + 5) 
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Indeed, e can be chosen as a truncation of the binary expansion of Note that 
without loss of generality we can assume that 

1 

since we can always enlarge n. Using these relations, we obtain 



\s{L,e)^t\ 



2" 



s(L) 



1 



,=0 



< t 



2"s(L) ' ^] 



1 _ 1 L Y^"-l ^ oj-n+l\ 

^ ''\2"s{L) ^ 2^3=0 ^J + ^^ ) 



2"s(L) 



< 



2"s(L) 
t2j 



1 



J. 



Note that for the last equality we used S < j. 

Now let e S Eoo be defined by Sj = for all j > jo for some jo E N, and let 
M > 0. Then there exists some n E N such that 

1 ""^ 1 
+ V Ei+i 2J-"+i < — for all n > uq, 

which implies 

(ei . . . £„)) = 1 , ,2^-»+i 

hence lim„^oo s{L, (ei . . . e„)) = oo. 

Finally, let e e i?oo be defined by Sj — 1 for all j > jo for some jo E N. Then, 
for all n e N, 

s{L, (ei . ..En)) = — 



and hence, 

lim s(i, (ei . ..£„)) = 



lim 

n — *C30 



1 

2"s{L) 



2-5^-5tEti (l-£,)2^- 



□ 



Thus only employing A/g and Mi we can move any line arbitrarily close to any 
line of slope € oo]. This shows the range of directions we might attain (compare 
Figure |3|. However, we would like to mention that the change of orientation of 
the data induced by the subdivision scheme (see Definition 3.2 1 is also affected by 
directionality of the masks. 

Theorem 2.7. Let L be a line through the origin with s E (0, cx)]. Then, for each 
t E [—-^, — oo] and S > 0, there exists some n E N and e E £"„ such that 

\s{L,e)-t\ < S. 



Similar results as Theorems 



2.6 



and 



2.7 



also hold for the matrices M,, e E 



{— 1, 0, 1}. We omit to also state these results for the sake of brevity, since they are 
similar to the previous theorems. 
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2.3. A Directional Refinement of the Lattice 1? . The results in the preced- 
ing section point out how to refine 1? in a directional way such that all possible 
directions can be attained. Dependent on whether we intend to map say the y-axis 
to a line with a slope contained in cx)], \—\, — oo], or [— 5, we choose to refine 
by using the matrices Mo,Mi, M_i,Mo, or M^i, Mq, Mi, respectively. Once the 
type of matrices is chosen, we iterate depending on the angle we would like to attain 



by using Theorem |2.6[ Theorem 2.7 or the corresponding result for the matrices 
Me, e e { — 1,0, 1}. For an illustration of the different areas of lines through the 
origin which can be attained during the refinement process dependent on the chosen 
matrices we refer to Figure |3] 




Figure 3. This figure shows the different areas of lines through 
the origin which can be attained during the refinement process 
depending on the choice of M^ and and ee{ — 1,0,1}. 



From now on we will focus entirely on the matrices Mg and Mi. All following 
results can be derived in a similar way for M_i, Mq and for M_i, Mq, Mi. 

3. Adaptive Directional Subdivision 

In this section, we finally arrive at the announced definition of a new type of sub- 
division schemes, based on the interaction of two "normal" stationary subdivision 
schemes, which we will study in the sequel. To that end, we choose two masks a^, 
e G {0, 1}, i.e., finitely supported sequences g ^oo (^^) well as the expanding 
scaling matrices = 1^1^^, e G {0, 1}. These matrices can be given explicitly as 

(3.1) 1^0= (0 fj and Wi^(^^ ' 

and again we set — Wg^ ■ ■ ■ W^^ , e £ £„ . Also note that 




Such a decomposition also exists for the iterated matrices W^, e G {0, 1}", 71 e N. 
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To formulate the next auxiliary result, we also define for e e i?„ the dyadic 
number 

n 

[£]2 = .£!...£„ :=^e,2-^' e [0,1]. 

With this notation at hand, we obtain the following counterpiece of Lemma 2.3 
Lemma 3.1. For n G No and e E En, we have 



where 



U.= \ 1 ^^= 1 



hence Ue — Vf 



Proof. The proof is again of inductive nature and relies on noting that 

A 0\ A" -4"2[e]2\ _ /4«+i -4"+i2 [(e,0)]2 
lo 2) I 2"+i j ~ I 2"+2 



as well as 



4 -4\ /4" -4"2[e]2\ /4"+^ -4"+^ (2 [e]2 + 2"") 

2 ) \ 2" ) ^ \ Q 2"+i 

4"+i -4"+i2 [(e,l)]2 
2"+i 



Hence, 

Since for a; e 1 

1 — x\ /I — fca; 



^0 1 y 1 
also the final claim follows. □ 



Note that V(o,...,o)i ^(1,0,. ..,0)1 ^^'^ are unimodular matrices, i.e., they have 

an inverse in Z^^^. A particular role will be played by the two matrices 




which satisfy 
(3.2) 

Wi = UWo = WoV, i.e. Wi = U-^WiV and Wq^UWqV-^. 

The associated subdivision schemes are now defined as follows. The term adap- 
tive refers to the tree-like structure, which provides various branches for subdivision, 
whereas the term directional refers to the directional structure which comes from 
the shearing process contained in the dilation matrices W^, e E E. 
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Definition 3.2. Let e £oo(^^), £ G {0, 1} be two masks, that is, two finitely sup- 
ported sequences, and let W^, e € {0, 1} be defined as in (3.1 1. Then the associated 



adaptive directional subdivision scheme of order n is defined by 



£ e En, n e N, 



where, for rj e {0, 1}, 

Sr/C ■■= Sa^.w^c := ^ (• - W^, a) c{a), c S ioo (Z^) , 

Note that both the mask as well as the scaling matrix of these subdivision 
schemes depend on the index e. Moreover, we wish to remark that these schemes 
can clearly be computed in a tree-like fashion by setting 

SeC = %',.„)C = S,„S,, = J2 (■ - W^^P) Se'ciP), S' € E„^i. 

Adaptive directional subdivision schemes can be considered subdivision schemes 
of their own, however, with a different scaling matrix. This is easily seen by means 
of the following example: for a G we have 

= ^ a,, (• - W,,I3) «ei - We,j) c(7) 



c(7) 



E 



pel? 



An inductive application of this argument immediately gives the next result. 
Lemma 3.3. For e G En, the subdivision scheme acts as 

S^c{a) = ^ ae (a - W^P) c(/3), a € Z^, 

where the coefficient sequences are recursively defined as = a(e',e„) — Se^a^i. 

To get a better understanding of the geometry of adaptive directional subdivi- 
sion, we write ai as ai = oq (C/-) which is always possible since U is unimodular. 



It then follows from repeated applications of (3.2 1 that 



Sai,WtC = ai (■ - Wio) c{a) 

= E "0 ■ -UWiU-^Ua) c{a) 

aeJ? 
aeJ? 

= Y'S.o{U ■ ~UWoV~^a) c{U-h 



SHEARLET MULTIRESOLUTION ANALYSIS 



13 



= '^aQ{U--Woa)c{U~^a) 

= {Sa,.,WoC{U-'-)){U-). 

This identity can be rewritten in terms of dilation operators as 

Si = Du So Djj-i = Du Sq D^^, hence S'(i^..._i) = Du <5'(o,...,o) ^c/^ 

and enables us to implement the subdivision scheme 5*1 in terms of Sq and the 
shear operator Du. Moreover, it explains the geometry of the scheme 5i: first, 
a shearing by U^^ is applied to the data sequence, then the subdivision operator 
refines the data in the sheared direction with a higher resolution than the data in 
the non-sheared direction, so that the additional application of the shearing by U 
does not fully compensate the initial one. In summary, this process leads to limit 
functions which are sheared versions of the limit function of 5*0 and the amount of 
shearing is determined by when and how often Si is applied in the process. We 
remark that this geometry is very much in the spirit of the Continuous Shearlet 
Transform, which can be regarded as applying a shearing operator, an anisotropic 
2-D Wavelet Transform, and again a shearing operator [23] . 

4. Convergence 

In this section, we shall study convergence of the previously introduced adaptive 
directional subdivision schemes. To that end, we introduce the projection operators 
Pn ■ ^ En, n g N, which extract the initial segment of order n from a 

sequence: P„e = (ei, . . . ,e„). 

Definition 4.1. The adaptive directional subdivision scheme is said to be conver- 
gent in C (M^), if, for any e £ Eao, there exists a nonzero uniformly continuous 
function /e G C (K^) such that 




Note that this is equivalent to 

lim sup |/e {Wp^ a) - ap^^{a) \ = 0. 

Since any sequence c e ^ (Z^) can be trivially written as 

c = ^ c(a) - a) , 6{a) -.^ 6afi, 

and since the subdivision operator is linear, we immediately obtain the following 
convolution style representation of the limit function. 

Proposition 4.2. I] the adaptive directional subdivision scheme converges for some 
£ £ Eaa then the limit function takes the form 

fe*c = ^ c(a)/e (• - a) . 
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4.1. Basic Properties. This definition of convergence has an immediate conse- 
quence: If the adaptive directional subdivision scheme is a convergent one, then, 
in particular, uq and ai must define convergent adaptive directional subdivision 
schemes, which follows by simply choosing e = (0, 0, ... ) and e = (1,1,...), respec- 
tively. Consequently, they must both preserve constants. 

Lemma 4.3. // the adaptive directional subdivision scheme is convergent, then 

(4.1) ^ ae(a-H We/3) = 1, a G Z^, e e {0, 1}. 

/3GZ2 

An alternative but equivalent definition of convergence of a adaptive directional 
subdivision scheme can be given in terms of function spaces instead of sequence 
spaces by means of test functions. 

Definition 4.4. A function g £ C (K^) is called a test function, if it is compactly 
supported and its integer translates form a stable partition of unity, that is, 

(i) = 

(ii) there exist constants < A < B < oo such that for any c G ioo 

A \\c\\^ < ||g*c||^ <S||c||oo, 9*c:= ^ c{a) g {■ ~ a) . 

The most prominent examples for test functions are the tensor product B-Splines 
so that there even exist refinable test functions of arbitrary regularity. With the 
help of test functions, convergence can be described as follows. 

Theorem 4.5. The adaptive directional subdivision scheme converges if and only 
if for any e € Eao there exists a nonzero uniformly continuous function such that 

(4.2) lim II/, - (5 * Sp^J) {Wp„e)\\oo - 

n— *oo 

(i) for some test function g. 

(ii) for any test function g. 

Proof. For classical subdivision, this result is due to Dahmen and Micchelli [TS] and 
we just show how it can be extended in a straightforward way to adaptive directional 
subdivision. To that end, let g be any test function and recall that for any uniformly 
continuous function / and any expanding matrix M the "quasi-interpolant" 

g*<JMf ^ ^ f{Ma) g{--a), 

with the sampling operator aM ■— {f{M a) : a £ Z^), satisfies 
11/ - 5 * aM-i/ (M.)IU <CgLo (/, ||M-i||) , 

where 

w(/,5):=sup sup \f{x)-f{y)\, 

a;eR2 \\x-y\\^<S 
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denotes the modulus of continuity of /. Recall that lo (/, (5) ^ for 5 ^ as long 
as / is uniformly continuous. Now, we have that 

\\h~{g*Sp^,5) (W^P„e-)lloo 



< 



< CgUj{f,\\Wp^\\\)+B 

On the other hand, 



"■w-i fe - Sp^i^S 



(^a^-iJ,-Sp,^j') (Wp,^,-) 

< (||/, - (g * <Jw-^Je) + ll/e - (.9 * ^P„s'5) (W^P„e-)lloo) 

which verifies the equivalence. Since therefore convergence of the adaptive direc- 
tional subdivision scheme is equivalent to (4.2) holding for an arbitrary test func- 
tion, this property holds for one particular test function if and only if it holds for 
any test function. □ 

Theorem 4.6. // the adaptive directional subdivision scheme converges, then the 
limit functions /e, e G £-00 , satisfy the refinement equation 



(4.3) 



/e = X! O-etia) fe{We^ ' , ? (£2, £3, • ■ •) 



Proof. We define the transition operator 

Tef^ E aM)f{W,--a), feC{R^), £€{0,1} 

and note that, for c € (-oo, 

{TJ)*c = E r,/(.-a) c(a)= ^ a,{(3) c{a)f {W, ■ -W,a - f3) 



By iteration, we then find for e £ {0, 1}" that 
(/ * S,c) {W,-) = (/ * S,,^ ■ ■ ■ S,,c) {W,,^ -...-W,,-) 

= (T,„/ * 5,„_, • • • S-e.c) •) = ••• = (7^/ * c) 

where 

TJ = T,,---T,J, £€{0,1}". 

Since, for n E N, 

T,J^^T,J^*5^{f^*S,,5){W,,-) 

= [{fe~ (.9 * ^P„_,e^) (M^P„_,e)) * S,,5] {W,, ■) 

+ [{g * Sp^_,?5) (W^P„_,e) * S,,5\ [W,,-) 
= [{f? - {g * Sp^_,^5) {Wp^_,^)) * S,,5\ {W,, ■) + {g* Sp^,5) {Wp^,-) , 
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it follows that 

\\TeJe-M^ 

< II (/e - (g * 5p„_,r) (M/P„_,s)) * ^e,<5||^ + ll/e - (5 * Sp^J) (M^P^eOIL 

< ll^sJI ||/s-(5*^P„-,e)(M^P„_,e)L + ll/e-(5*5p„e<5)(T4^P„e-)IL 

and the right hand side of this inequality converges to zero for n — > cxd while the 



left hand side is independent of n. Thus T^^^fe ~ fe which is (4.3). □ 
4.2. An Algebraic Description, Sum Rules and Polynomial Reproduction. 



Next, we give a more detailed description of the necessary condition (4.1) from 
Lemma |4.3| in algebraic terms. To that end, we recall the definition of the symbol 
of a mask a, defined as 



as well as the subsymhols 



zeC2 = (C\{0})' 



T^e He-.^ wj [0,1)^ nz^ 



{0,1}. 



The symbol can be "reconstructed" from the subsymbols by the well-known formula 



a* (z) = 5^ z^' a* 
from which the following result follows immediately, cf. 



ee {0,1}, 



Proposition 4.7. The mask satisfies (4-.1), the sum rule of order 0, if and only 

a*{z) = 0, z e |e-2-^^-"^" : 1] ^ H,\ {0}} . 

For a more algebraic description, we need the notion of a quotient ideal. Recall 
that an ideal in A, the ring of Laurent polynomials in two variables, is a subset of A 
that is closed under addition and multiplication by arbitrary Laurent polynomials. 
The quotient ideal of two Laurent ideals /, J, is defined as 

/: J:={/e A : / • J C /} 

and has the almost obvious property that I C I : J. For any matrix X e Z^^^, 
with column vectors xi,X2 we finally define the ideal 

(^''-l) :=(z^^-l,z^^-l) :={/i(z)(2^^-l) + /2(z)(z-^-l) : /i,/2eA} 

and its special case (z — 1) := — l). Then we have the following result from 



Theorem 4.8. The mask satisfies ihe sum rule of order 0, ij and only if 

a* e (z^- - 1) : (z - 1) . 

To conveniently formulate an important consequence of this theorem, we introduce 
the vectors 



[z- - 1] 



Z^l - 1 

^X2 _ I 



72x2 



With this notation we have the following result. 
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Corollary 4.9. // the adaptive directional subdivision scheme converges, then there 
exist matrix valued masks B^, e G {0, 1} such that 

(4.4) [z - 1] a^iz) = B;{z) - l] , e E {0, 1}. 

Proof. Any convergent subdivision must satisfy the sum rule of order for a^, 



£ e {0, 1}, and so, by Theorem 4.8 it follows for e e {0, 1} and j — \,2 that 

(z, - 1) <(z) = b*,{z) (z(^-)i - l) + b*,{z) (z(^^)= - l) . 
Written in matrix form, this is what has been claimed. □ 



Definition 4.10. The matrix masks B^, e e {0, 1}, from (4.4) are called represen- 
tation masks of a^, £ G {0, 1}, respectively. 

Remark 4.11. Recall that the computation of the representation masks can be 
performed by reduction, a multivariate generalization of division with remainder, see 
[To. 2lJ ^ for the term order and homogeneous versions of this process, respectively. 
Therefore, the symbolic determination of B^ can easily be done with the help of 
practically any Computer Algebra system that supports constructive polynomial 
ideal theory. 

Note however, that the representation masks are not unique to the appearance 
of syzygies of — l] , not even if an H-representation, cf. [25J, is chosen where 
- in the case of Wq - we have the "minimal degree" requirements that 

deg 6ii = deg 621 = deg oq - 3, deg &12 = deg 622 = deg oq - 1, 

see also l28l. 



We continue by giving explicit bases of the quotient ideals for our specific choice of 
We. This is easy for Wq as all entries in this matrix are nonnegative, and indeed it 
is not difficult to see that 

lo := (z^«-l):(z-l) = (zf-l,Z2'-l):(^-l> 

= {{zl + zf + Z, + l) (Z2 + 1)) + {zt - 1) + (z| - 1) . 

In fact, the graded homogeneous leading terms of the above ideal basis are zlz2, z\ 
and z| so that the quotient space is spanned exactly by the seven monomials 

I, Zl, Z2, zl, Z1Z2, zl, zlz2, 

and their number coincides with the number of joint zeros of /q. Hence, by the 
same reasoning as in [551 [M| they even form a graded Grobner basis, hence an 
H-basis of the ideal Jo . Recall that a subset H of an ideal / is called an H-basis, 
if any polynomial f (z I can be written in the form 

f =^ fhh, deg / > deg + deg h, 

heH 

where deg denotes, as usual, the total degree of a polynomial. We will also use n„ 
for the vector space of all polynomials of total degree at most n. 

The situation for Ii — (^z^^ — l) appears to be a little bit more intricate due 
to the appearance of a negative entry in W\. Here it is helpful to recall that 

W\ — L/Wo, ^ = ( Q ^° define y = z^ = (zi, Zj"^Z2), hence also z = ^ — 
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(2/1 1 2/12/2) and to realize that 

h = (z^^-l):(z-l>^(z^^»-l):(z-l> = (y^"-l):(/^"^-l) 
Since 

(2/^ ' - 1) = (2/1 - 1,2/1272 - 1) = (2/1 - l,2/i2/2 - (2/12/2 + 2/2) (2/1 - 1) - 1) 

= (2/i-i,2y2-i) = (2/-i>, 

we thus obtain that 

h = (2/^" - 1) : (2/ - 1> 

= ( (2/? + 2/? + 2/1 + 1) (2/2 + l)) + {yt-l) + {vl - 1) 
= {{zl + zl + Zl + 1) (z2 + zl)) + (zf - 1) + (z2 - zt) 

To arrive at the somewhat surprising observation that in fact Ii — Iq, we add 
zf — 1 to the third basis element, z| — zf, yielding z| — 1 again, and subtract 
(zi + 1) (zf — 1) from the first basis element which leads to 

{zf + zf + Zl + 1) (Z2 + zf) - (zi + 1) (zf - 1) 

= {zf + zf + zi + l) Z2 + zf + zf + zf + zf- zf- zf + Zl + 1 

= {zf + Z? + Zl + 1) (Z2 + 1) 

and therefore to the following result. 

Theorem 4.12. The two quotient ideals = (^z^^ — l) : (z — 1), e G {0,1}; 
coincide and have the H-basis representation 

(4.5) / := /o = /i = {{zf + zf + Zl + 1) (Z2 + 1)) + {zf - l) + (z^ - l) . 

The fact that Iq — Ii may appear a little bit surprising at first view, since it implies 
that, for any finitely supported mask a, we have 

^ a(a+ VKo/3) = 1, a e ^ ^ a (a + T^i/?) = 1, aeZ^. 

Hence the necessary "sum rule" condition with respect to Wq is equivalent to the 
one with respect to Wi. However, if we write Wi = WqV with the unimodular 

matrix V = ( ^ 1 M , then a simple change of the summation variable indeed 



^0 1 
gives for any a G 



^ a (a + WiP) = ^ a (a + WqVP) = ^ a (a + Wof3) 

;3eZ2 /3eZ2 



and confirms ( |4.5[ ). 

Moreover, note that Theorem |4.12| gives a way to parameterize the ideal of all 
admissible polynomial masks. Indeed, for any n G N we have that 

/ n Hd = p{z) {zf - 1) + q{z) {zf + z2 + z + 1) (Z2 + 1) + r(z) (z| - l) , 
with degp < n — 4, deg q < n — A, and deg r < n — 2. 

For a polynomial of this form, the decomposition with respect to Wq, i.e., the 

matrix polynomial Bq, becomes 

(4.6) 

. /(zi-l)p(2) + (z2 + l) <z(z) (zi-l)r(z) \ 

" ^ ' ~[ {Z2-1) p{z) {zf + z? + Zl + 1) g(z) + (Z2 - 1) r(z) .) 
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Since the two Laurent ideals Iq and Ii coincide, the decomposition of a^f into 
takes exactly the same form as B'q in (4.6 1. 



Next, we rephrase the identity (4.4 1 by means of the backwards difference operator 
V, defined for a sequence a as 

^-=(:(::S)::(:))' (v«rw = [-i]«*w- 

where 771 — (Y\ and 772 ~ (^] denote the unit multiindices in Z^. Since, in 



addition, any finitely supported matrix sequence B satisfies 

{Sb,w^c)* (z) = B*{z) c (z^=) , c e V^oo(Z2), e G {0, 1}, 

where 

Sb,w^c:= J2 Bi--W,a) c(a), 

our quotient ideal representation ( |4.4[ ) can equivalently be written in terms of the 
difference operator as 

V5a,,M/. = 5b.,W,V, £€{0,1}. 

We end this section by recalling that quotient ideal containment also charac- 
terizes the order of polynomial reproduction provided by the two masks and thus 
the subdivision scheme. Recall that a mask a provides polynomial reproduction 
of order n, if the leading forms of all polynomial sequences are reproduced by the 
scheme: 

Sallk = life, fc = 0, . . . , n, life := < ^ a'^ : a e Z^ 

[h\<k 

Polynomial reproduction is essential for the smoothness of the refinable limit func- 
tion fB] as well as for the approximation order of the associated wavelet construction. 
With the methods from |26l I28j we can now easily describe polynomial reproduc- 
tion. 

Theorem 4.13. The directional subdivision scheme preserves polynomials of degree 
n, i.e., S'sIIfe = Ilfc, e <E E, k — 0, . . . ,n, if and only if 

4.3. A Characterization of Convergence. Finally, we will give a characteriza- 
tion of convergence of the adaptive directional subdivision scheme, like usually in 
terms of a (restricted joint) spectral radius. In this subsection, the adaptive direc- 
tional subdivision scheme both for masks ag and ai as well as for their associated 
matrix sequences Bq and Bi will come into play. To distinguish both, for the first, 
we again employ the notation S^, s £ E, whereas the second adaptive directional 
subdivision scheme will be denoted by , e £ E. 

Now, given two matrix masks B^, e £ {0, 1}, their restricted joint spectral radius 
is defined as 

p (i?o, -Si I V) = limsup sup sup ||'S'^c|| 

ri^oo ee{oa}" ceV^oo 

The joint spectral radius is called "restricted" since the supremum is not taken over 
all 2-vector valued sequences but only over the proper subset V^oo , see [HI [30] . The 
main result of this paragraph is now as follows. 
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Theorem 4.14. The adaptive directional subdivision scheme based on the masks 
flg, e G {0, 1} converges if and only if a*(z) £ I and the representation masks B^, 
ee {0,1}, satisfy p {Bo, Bi\V) < 1. 



We will split the lengthy proof of Theorem 4.14 into several partial results, begin- 
ning with the sufficiency of the spectral radius condition. To that end, we will show 
that, starting with a particular test function g, the sequence g * Sp^^^c converges to 
a limit function for any choice of e € i?oo and any c G ioc- Indeed, we choose the 
test function g to be Wo-refinable with respect to a mask b, that is 



(4.7) 



9= 3 {Wo ■ -a) 



Such functions can be easily shown to exist, even with an arbitrary order of smooth- 
ness: pick any cardinal B-spline cj) — M {■ \ 0, . . . , N) with refinement mask h, then 
a double application of the refinement equation with respect to the first variable 
shows that the tensor product function 

9{x, y) = [{(j) 0] {x, y) = (0 * (f) (x) (t){y) -. ip{x) (j){y) 

is Wo~refinable with respect to the mask b = Sh,2h ® h, where Sh,2 denotes the 
subdivision scheme with mask h and dilation 2. The following lemma states a more 
general process. 

Lemma 4.15. Let 6i, 62 G ^(Z) be 2-refinable masks, and let the mask bi be defined 
by 61 (to) = S'f,i,2&i(TO) ~ X^feez ^1(^)^1 ("^ ~ 2^)- Then the mask cq = 61 (8) &2 is 
Wo-refinable, and ai = ao{U-) is Wi-refinable. 

Proof. Let Lpi, ip2 be univariate functions which are 2-refinable with respect to b\,bi, 
respectively, i.e., 

^^ = Yb,{k)ip{2■-k), i = l,2. 
fcez 

We claim that the function / defined by 

/ = V3l (g) (/32 

is Wo-refinable with respect to oq. Indeed, for x — (xi, 0:2) € K^, we obtain 
{bx®b2){a)f{Wox-a) 



E 




61 (ai - 2fc)6i(fc) ^i(4xi - ai) 



fcez 



2^6i(fc) 2^ &i(ai)(/7i(4xi 



E 

.Q2ez 



b2[0L2)V2{^X2 - 02) 



2fc-ai) 



•^2(X2) 



^ 61 (fc) (^1(2x1 

Lfcez 



fc) 



^2{X2) 



The claim concerning Wi-refinability of a\ follows from Lemma 4.16 



□ 



There also exists a canonical Wi-refinable function associated to g. 
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Lemma 4.16. If go — 9 is Wo-refinable with respect to the mask bo — b, then 
9i — 9a {U-) is Wi-refinable with respect to the mask bi = bo (U-). 

Proof. Setting gi ~ go {U-) and thus go = gi (U^^-), we find for a; G M'^ that 
9i{x) = 9o{Ux)^ J2^o{a)9o{WoUx-a)^ J2^o{a)goiWoV^x-a) 

= ^o{a) gi {U-^WiVx - U-^a) = ^ &o {Ua) gi {Wix - a) 

= ^ 6i(a) gx {Wix - a) , 

hence gi is Wi-refinable with respect to 6i. □ 

The next two observation are again of a more algebraic nature. 

Lemma 4.17. Suppose that a mask a satisfies Sa.w^c — for all constant se- 
quences c and some e G {0, 1}. Then there exists a 1 x 2 matrix mask B such that 

Proof. Again we refer to [261 128] where it has been shown that Sa.w^c — for all 
constant sequences c if and only if a*(z) e (^z^^ — l) which is in turn equivalent to 
the existence of a representation 

a*{z) = bl{z) (z('^-)i - l) + (zf'^-)^ " = ^* W \.''^' - 1] ' 

which is nothing but Sa,w^ = Sb.w^'^ ■ D 

Lemma 4.18. Suppose that a compactly supported function f satisfies f * c = Q 
for all constant sequences, then there exists a compactly supported, continuous 1x2 
matrix function G such that f * c ^ G * Vc for all c £ £oo (^^) • 

Proof. For any x G [0,1]'^ we consider the sequence fx = (/(x + a) : asZ^). 
Since / is compactly supported, any such sequence fx, x G [0, 1]^ has finite support 
and since / is continuous, the map a; i— > /^^ is a continuous one. 

By assumption, /a; *c = for any x and any constant sequence c, hence, with the 
scaling matrix /, the same methods as above yield that /* G (z^ — l) = (z — 1). 
Consequently, we have that 

/:(^) = 9:.l{z) {Zl - 1) + gl2{z) {Z2 - 1) = G;(z) [z - 1] 

where, like f^ and /^(z), also G^iz) depend continuously on x as they can be 
obtained by applying the orthogonal reduction process from |27j . Therefore, the 
function G, defined as 

G{x + a) = Gx{a), x G [0, l]^ a G 

has the properties claimed in the statement of the lemma. □ 

Now we are in position to prove the sufficiency of the spectral radius condition 
which we state as a separate proposition. 

Proposition 4.19. The adaptive directional subdivision scheme based on the masks 
Qe, e G {0,1} converges, if a*^{z) G / and the representation masks B^, e G {0,1}, 
satisfy p {Bo, Bi\\/) < 1. 
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Proof. For any 9 G (0,1 — p), there exists, by standard properties of the (joint) 
spectral radius, a constant C > such that 

||Sf„,Vc||^<C(p+0r = Ca", neN, e e E^, c ^ {Z") , 

where Q < a := p + 8 < 1. 

Now, let e e Eao be given and suppose first that £„ — 0. Then, by the refinability 
of the test function g from (4.7 1 and Lemma 4.17 which ensures the existence of a 
finitely supported matrix mask F such that 5o — Sb.Wa ^ SF,Wa^ , have that 

\\g * Sp^^c (W^p„e-) -9* Sp^_,^c {Wp^_,e-) 11^ 

= h * Sp^eC{Wp^e-) - g* Sb,WoSp„_isC{Wp^e-)\\^ 

= llg * {So ~ Sb,Wo) Sp^_iec\\^ 



< Bg \\{Sf) — Sb^Wo) Sp„_iec\ 



oo- ^9 \\SFMo^Sp^_^eC\\ 



SF,WaSp^_^e 



Vc 



< Bn 



If on the other hand £„ = 1, by using the function gi — giU ■) (cf. Lemma 4.16) 
we pass to the estimate 

\\g * Sp^eC (W"p„e-) - 3 * 5p„_,£c (T4^p„_ie-) 11^ 
< IKff - 9i) * 5p„ec(M^P„e-)IL + 11(5 - 9i) * 5p„_,,c(M^p„_,,.)L 

+ ||.9l * Sp^eCiWp,^e-) - 51 * '5'p„_iec(W^P„_ie-)lloo ' 

For the first two terms we now make use of Lemma [4. 181 to obtain that 

ll(5-5l)*^P„eC(M/p„e.)lloo 



= 11(5 - 9i) * ^P„.c|U = l|G * V5p„,c|U = \\G* 5|„eVc|L < Ca" 



and 



\\{9-9i)*Sp^_,,c{Wp,^_,,-)\\^<BGCa--\ 

respectively, while the third term can now be estimated as above again. In summary, 
we obtain that there exists a constant D > such that 

\\g * Sp^eciWp^e-) - 9 * Sp^_,,c{Wp„_,,-)\\^ < D a''-' 

so that for m G N 



19 * Sp^ 



■{Wp^ 



9*Sp„,c{Wp„,-)\\^<D 



1 



In other words, the sequence g * S'p^^c (VFp^g-) is a Cauchy sequence of continuous 
functions and thus must converge to a limit function for n oo. Convergence of 
the subdivision scheme then follows by standard means. □ 



The proof of the converse statement of Proposition 4.19 is based on the estimate 



- l|V5p„e<5|| 



Sp„eS (• - m) - Sp„eS (•) 
Sp,,eS (• - V2) - Sp^J (•) 



max \\Sp^,S (• - 7]j) - Sp^.S (OH, 



< max( S'p„e<5(--?7j)-cr^-i /(--Tyj) + Sp^J ~ cr^^-i f 
'^w-\fi--Vi)-'^w-'fi-) )' 
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hence, 

<max(2 Sp^.S - (Jy^-i J + <y^-i f {■ - rij) - a^-i f {■) ). 

If we assume that the subdivision scheme converges with uniformly continuous hmit 
function, then the right hand side converges to zero, hence also \\Sp ^ Vc||^ — > 
for n OO and any c € £00 (^^)- This, however, is not sufficient for our purposes. 
To show that the restricted spectral radius of p{Bo,Bi | V) is less than one, we 
have to show that 



(4.8) 



5i.Vc|L <C0„ II Veil 



lim Or, = 0, 



which will be prepared in the next lemmas. Here we follow the outline of a proof 
from [6j and show that there exists a constant C > such that 

||V5p„ec||^ < Cf?„ ||Vc||^ , lim 0n = 0, 



from which (4.8 1 follows immediately. We begin with an estimate on the limit 
function f^. 

Lemma 4.20. // oq and ai define a convergent subdivision scheme, then there 
exists a constant Ci > such that for any e G £^00 and any c e £00 (2^) 

\f,*c{x)- f,*c{y)\<C\Luif„S) \\Vc\\^, llx-yll <^ <1. 



Proof. Since, according to Lemma |4.3| convergence implies the preservation of con- 
stant sequences by the subdivision scheme, we also have that 



and thus, for any c € £00, any w G M and any x, y € with ||a; — y|| < S, 



\fe * c{x) - fe * c{y)\ = 



{f,{x - a) - fe{y ~ a)) (c(a) - w) 



where 



< #f^a;.i/ • uj{fe,^) ' max |c(a)-w| 



n^,y^{aeZ^ : f,{x-a)^0}u{aeZ^ : fe{y-a)^0} 



Since is finitely supported, we have that ^^x.y < 00. Specifically, if we assume 
that /e is supported on [-iV, iV]^, then i^fl^^y < (27V + 2)^ as long as S < 1. 
Choosing 

w = - { max c(a) + min c(a) 



2 \aen, 

it follows for any a € ^x,y that 
1 



\c{a) ~ w\ < 



max c{a) + min c{a) 



< II Veil 



hence, 



|/e*c(a;)-/,*c(y)|< i(27V + 2)^u;(/„<5) ||Vc||, 



as claimed. 



□ 
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The next result concerns the difference between the subdivision scheme and the 
limit function. 

Lemma 4.21. // the adaptive directional subdivision scheme based on the masks 
tte, e € {0,1}, then there exists a constant C2 > such that, for any n G N, we 
have 

Proof. We fix n, set, for abbreviation, 'e — Pn£, and assume again that as well 
as oo and ai are supported on [—N,NY. Again, we make use of the fact that 
and preserve constant data and obtain, for any a e and w G M, that 

S^c{a) -f,*c (Wf'a) = ^ {a^{a - Wr/3) - f, (Wf'a - /3)) (c(/3) - w) . 
Since 

f^a,?={aeZ2 : {a - Wef3) ^0}U {a eZ^ : fe (Wf^a - P) ^ 0} 

again satisfies ^^a.e ^ (2-/V + 2)^, the same judicious choice of w as above leads 
to the estimate 

\S^c{a)~ f,*c{WZ^a) \ < (27V + 2)* sup \a^{a) ~ f, {W^^ a)\ ||Vc||^ , 

from which the claim follows immediately. □ 

Now it is easy to complete the proof of the converse statement for convergence 
which we formulate in the following way. 

Proposition 4.22. // the adaptive directional subdivision scheme based on the 
masks a^, e G {0,1} converges then a*{z) G / and the representation masks B^, 
£G {0,1}, satisfy p {Bo, B^\V) < 1. 



Proof. In Lemma 4.3 it has already been shown that convergence implies a*(z) G /. 
Moreover, Lemma 4.20 and Lemma 4.21 allow us to conclude with C — max{Ci , C2} 
that, for any c G ioc, we have 

< II V (5p„c- /, * c(T4^p^V))L + II * ^(^p„V)L 

< 2c{\\Sp^J-f4Wp^\-)\\^ + u;{f,,\\Wp^\\\)) \\Vc\\^, 

and since 

lim \\Sp^J - f, {Wp\-)\\ + CO {f„\\Wp\\\) = 

by convergence of the adaptive directional subdivision scheme and uniform conti- 
nuity of the limit function, our prove is complete. □ 

5. Numerical Experiments 

In this section we present some numerical experiments which illustrate the ability 
of the developed class of subdivision schemes to adaptively change the orientation 
of the data. 

First, we recall that there exist a general way to construct masks, which are 
refinable with respect to the dilation matrices Wq and Wi, compare Lemma [4.15[ 
Now let the mask b G ^(Z) be chosen by 5(-3) = -3^ = 6(3), 6(-l) = ^ = 
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6(1), 6(0) = 1 and 6(771) = otherwise, which coincides with the mask studied by 
Deslauriers and Dubuc [15 . We remark that this mask yields a 2-interpolatory 
subdivision scheme (compare also Section [6]). By Lemma 
60 6 is Wo-refinable, and ai — ao(J7-) is Wi-refinable. 
In Figure |4] we illustrate the refinement of the matrix 



4.15 we know that oq 



(5.1) 



Ci = 







(a) 



1000 2000 3000 4000 5000 6000 7000 8000 9000 



(c) 



(d) 



Figure 4. This figure shows the refinement of the matrix Ci 
defined in after applying S, with (a) e = (0,0,0,0,0), (b) 

e = (0, 0,0,1,0), (c) £=(0,1,0, 0, 0), and (d) e = (0, 1, 1, 1, 1). 



and in Figure [5] we subdivide the data given by 



(5.2) 



C2 



/O 



1 







1 



0/ 



In both figures we employ different iterations of the subdivision schemes Sq and 
^i. As can clearly be seen, the application of 5*1 increases the angle the resulting 
images is sheared in the x-direction, where the angle depends on the particular 
path in the binary tree (see Figure [2]) we choose. 
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1000 2000 3000 4000 5000 6000 7000 SOOO 9000 10000 



1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 



(a) 



2000 4000 6000 3000 10000 12000 



(c) 



(b) 



4- 



2000 4000 



aooo 10000 12000 



(d) 



Figure 5. This figure shows the refinement of the matrix C2 
defined in ([K2| after applying with (a) e = (0,0,0,0,0), (b) 
e = (0, 0, 0, 1, 0), (c) e = (0, 1, 0, 0, 0), and (d) e = (0, 1, 1, 1, 1). 



6. Shearlet Multiresolution Analysis 

In this section we will show how the adaptive directional subdivision schemes de- 
veloped in the previous sections can be applied to derive a shearlet multiresolution 
analysis. For the sake of simplicity, in the computation of "dual functions" we will 
restrict ourselves to interpolatory subdivision schemes in this paper. Our idea is 
inspired by similar ideas for the construction of a fast wavelet decomposition from 
interpolatory subdivision schemes [17 . The construction of a shearlet multireso- 
lution analysis associated with general adaptive directional subdivision schemes is 
beyond the scope of this paper, and will be studied in a forthcoming paper. 

Before constructing the scaling spaces we first need to discuss whether there exist 
masks ao and a\ such that the subdivision schemes 5*0 and S\ are both interpolatory, 
respectively, which immediately implies that is interpolatory for each e £ -Eoo- 
To that end, we proceed by using a tensor product approach. Recall that a mask 
ao leads to an interpolatory subdivision scheme Sq provided that 

(6.1) aoiWoa) = 6^,0 for all a e 1? , 

likewise does a mask a\ lead to an interpolatory subdivision scheme S\ provided 
that 

(6.2) ax{Wxa) = b^s) for all a £l? . 



There exists a canonical way to define Oi by means of the matrix U as indicated 
by the following lemma (compare also Lemma 4.151. 



SHEARLET MULTIRESOLUTION ANALYSIS 



27 



Lemma 6.1. Let foi, 62 G ^(Z) be masks which satisfy bi{2m) ~ 6^.0 for all m G Z, 
i — 1,2 and let the mask hi be defined by 61 (to) = Sbi,2bi{m) — '^k^j^bi{k)bi{m — 



2k). Then the mask hi (E) 62 satisfies (6.1 1, and the mask (hi (g) b2){U •) satisfies 
K2\. 



Proof. Given some a = (ai,a2) G Z^, we obtain 
{hi (g) b2){Woa) = ^ 6i(fc)6i(4ai - 2k) ■ 620,2,0 = ^ bi{k)S2a,-kfi ■ Sa2,o = ^a,o- 

A similar computation shows (5® b){UWia) = Sa,o. □ 

Suppose we have chosen masks oq and ai so that the subdivision scheme is 
interpolatory and converges for each e G E^o ■ To define the scaling functions, recall 
that we wrote e* = (£,0,0...) for the canonical embedding of E into E^o] the 
image of this embedding operation, 

E* ^{e* : e e £;} C E^ 

thus consists of all infinite 0-1-sequences which contain only a finite number of 
nonzero components. It is worthwhile to keep in mind that the subdivision scheme 
converges for all £ G -E* if and only if ag defines a convergent subdivision scheme 
and hence the functions 

{/e : £ e E*} = {/,. : sgE} 

which will be needed to build the MRA can be ensured to exist by requiring the 
existence of an appropriate solution of the refinement equation associated to ag. 
This is a much weaker condition, of course, than convergence of the for any 

£ S Eao- 

Definition 6.2. The shearlet scaling spaces are defined as 

Vb = span {/e. (• - a) : aGl?, e G E] 

and 

K = '^e.n> 1, 

eG{0,l}" 

where 

Ve = span {/ (VFe --a) : aeZ^, f e Vq} for all e e E. 

Indeed this choice of scaling spaces provides a multiresolution analysis, which 
is the focus of the following theorem. The main ingredient in the proof is - as it 



should be - the refinement equation (4.3 1 



Theorem 6.3. The spaces (Vn)n>o create a multiresolution analysis. In particular, 

(i) the spaces Vn, n >0 are translation invariant, 

(ii) Vn ^ Ki+i for all n > 0, and 

(iii) for each n e N, we have f € Vn <^ fiW^ ■) e Vn+i for each e G {0, 1}. 

Proof. Statement (i) follows immediately from the definition of Vn, which is a trans- 
lational completion. 

To verify the nestedness property (ii), we consider an arbitrary "basis element" 
f EVn of the form 

(6.3) f = fr {W, ■ ~a) , eeEn, tj ^ (rj^,^) e E, a G Z^, 
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and make use of the refinement equation ( |4.3[ ) to verify that 

/ = ^ a^,(/3) ff}, (W^.iW, ■-a)-f3)=Y, " fv' i^,, ■ -/3) 

with e' = {e, rji) e En+i, hence f E Ve' C Vn+i- 



To verify (iii) we again consider a function element f £ Vn of the form (6.3 1. 
One implication follows from 

f{Wr-)^fr{W^r,e) - -a), re{0,l}, 

the other one can be deduced in a similar way by considering / £ Vn+i and showing 
that this yields / [W''^- ) G F„ for any r e {0, 1}. □ 

Notice that for each fixed e £ E, the set of functions /e* (• — a), a G Z^, can be 
interpreted as being derived from Sa by refining with the subdivision scheme S^- 
Since is interpolatory, this set of functions is linearly independent. 

Some of the scaUng functions which generate Vq are plotted in Figure [6] The 
different orientations due to the application of the adaptive directional subdivision 
scheme to the Dirac delta i5o is evident. This fact forces the associated shear let 
spaces to also comprise directionality, hence to react to directional behavior of the 
data. 




(c) (d) 

Figure 6. This figure shows the refinement of i5o after applying 
with (a) e = (0, 0, 0, 0, 0), (b) e = (0, 0, 0, 1, 0), (c) e = (0, 1, 0, 0, 0), 
and (d) e = (0,1,1,1,1). 
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7. Fast Shearlet Decomposition 

Let Vm n G No, denote a sequence of projections from Vn+i to Vn, respectively, 
and define the shearlet spaces as = {Vn ~ I) Vn+i, n G Nq, hence as an appro- 
priate complement of Vn in Vn+i- In classical MRA, V is chosen as an orthogonal 
projection, but following the approach from fT8], we can also use interpolation as 
a projection, provided that the subdivision schemes were interpolatory. 

7.1. Refinable Functions. In order to establish the shearlet decomposition, we 
require the following two observations. 

Lemma 7.1. For all e £ E and c £ ^(Z^), we have 

^ c (W-'a) fo (W, . -a) = ^ c (W.-^a) fo {U, {W^ ■ -a)) . 

Proof. Since all the matrices U^, e £ E, are unimodular, we obtain 

^ c (W-'a) fo {W, --a) = ^ c (Wo^'U^'a) fo {W, ■ -a) 

= ^ c(W^o-"a) /o(C/eW ■-«))■ □ 

To formulate the next result, we denote r : E —* E the reversal operator for 
sequences, which maps e — (£i,...,£„) to r{e) := r(£i,...,£„) := (£„,..., £i). 
Moreover, we will write 0^ = ^a;0* for the zero sequence in E^, k £ N. We can now 
derive the following crucial relationship between refinable functions and subdivision. 

Lemma 7.2. For < k < n, e — {ri,T) £ E, rj £ E^. and c £ we have 

J2 c{a) /e. (Vt^o""' ■ -«) = E SA^)fr' {Wriv)W--' ■ -a) . 

Proof. Without loss of generality we can assume that t — (0). Then, for £ = (£i, £), 
the refinement equation ( [4.3^ gives 

This is the initial step for the inductive proof that for j < k we have 
(7.1) E^H/-*(^o""'--") 
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Indeed, applying the refinement equation (4.3) once more to (7.11, we get that 



which advances the induction hypothesis in (7.1 1. Specifically, for j = k this identity 
gives 

^ c{a) (l^o"-''^ • -a) = J2 ^-^W /o (W^K.,o„-.) ' 

= ^ S,c{f3) h {UrieM,._,)W^ ■ -(3) 

= J2 (t^'-(e,o„-.)/3) /o (C^.(e,o„_.) (W^o" • -/?)) ■ 

Since for any rj e 

fe 

-2"+Mr(77,0„„fc)]2 = -2"+i2-"+'=^77fe_,2-^' 

i=i 

fc 

= -2'=+i^r(r;),2-^- 

= -2'=+MK'?)]2> 

we finally get the identity 

which proves the claim. □ 

Now suppose we are given some data from a finely sampled function on the grid 
VFg~"Z — 4~"Z X 2~"Z, say. The key idea for the decomposition of this data, 
dependent on different directions, is stated in the following result which is the 
backbone of the MRA based fast discrete shearlet decomposition. We would like 
to mention that it relies on the fact that the masks ao and ai are chosen to be 
interpolatory and thus give us an explicit expression for Vn ~ I- 

The wavelet part of such a decomposition is, as usual, related to the represen- 
tatives of the quotient groups :— Z^/Wr(e)Z^, e € E. Since for e G En, we 
have det Wq = det Wi = 8" , all such quotient groups consist of a number of ele- 
ments that depends only on the length of e; we will denote by V* a selection of 
8" — 1 representatives for \ {[0]}. In the sequel, we will make use of the notation 
Dmc — c(M-), M being some 2x2-matrix. 
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Theorem 7.3. For c e t{I?), e ^ {ri,T) e E, rj e Ek and n > k we have that 

(7.2) + ^ ^ (c - 5,i^vy,(,,c) + 7) • -a) - 7) . 

Proof. The decomposition is based on the prediction-correction method which has 
become standard for interpolation based wavelet decomposition, in particular in 
connection with the so-called "lazy wavelet" and the associated "lifting schemes" 

m- 

We subsample the data c e £ (Z^) to obtain c' — Dw^^^-^c and make use of 
Lemma [7.21 to obtain that 



This identity is then decomposed with respect to giving the prediction 

+ E E + 7) /r* (W^.C)) (W^O""' • -«) - 7) 

7er; Qez2 

since the subdivision schemes were supposed to be interpolatory. Comparing this 
with the decomposition 

E C{a) fr. {Wr^rj)W^-'' ' "«) = E E (WrM^o'' ' " l) 

= E^(")/--(^-w(^o""'--")) 
+ E E ^(") i^ri,) {wr" ■ -a) ~ 7) 

7er* aGZ2 



we have to apply precisely the correction from (7.2 1. □ 



For the special case rj = ei and thus t = e , Theorem 7.3 simplifies into the following 
form. 
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Corollary 7.4. For c <E £ {I?), e e E andneN we have that 
c{a) {W,,W--' ■ -«) = E ^(^-i") (^o""' ■ 
(7-3) + E E - Se.^w., c) (M/,,a + 7) /e> (W^e, (W^o""' • -a) - 7) 



Remark 7.5. The decomposition (7.3 1 is the shearlet decomposition associated with 
the shearlet MRA: The function on the left hand side belongs to Vn and is written 
as the sum of a function in Ki-i and correction terms from Vn that vanish at Wi^-^lP' 
— the shearlets in the interpolatory MRA. 

7.2. Decomposition Algorithm. The fast shearlet decomposition is now based 
on an iterative application of (7.3 1, where each step can be understood as filtering 
by means of a filter bank. To that end, we have to interpret the initial sequence 
c G (Z^) appropriately Denoting by := /o(C^e-) the "sheared" version of the 
refinablc function /o, we form the quasi- intcrpolants 
(7.4) 



These are precisely the functions which appear on the left hand side of (7.2 1 and 
(7.3 1. It is worthwhile to note that all the functions Qe^n are relying on the same 
initial data c e £ (Z^). 

The interpretation of ( 7.4 1 is rather easy now if we take into account that /o 
was assumed to be the limit function of an interpolatory scheme, hence cardinal: 
fo (a) — (5o,Q, a G Z^. Hence, since 

(7.5) ge,«(a;) = E ■^o (^-^ - ' xeR\ 

we can substitute x = W~^a = M^a and use the cardinality of fo to find that 
(7e,„ (MgOf) = c{a) or q^^n (Wq^o) = c{Uea), respectively. The latter tells us 
that we should interpret the sequence c as a function sampled at the grid Wq"^!? , 
while the parameter e determines how this data is sheared and which thus arc the 
directions "preferred" by the wavelet decomposition. 

For the fast decomposition we now start with c £ £ interpret it as in (7.51, 
and decompose it in two ways, namely, for e G into 

qe.n = E (") (^o""' ■ -«) + E E ^^-.7 («) /O [^e (t^o""' ' ~ l) : 

where the coefficients 

Ce = Dw,c 
4,^ = {c-S,Dw^c){W,-+-f) 

are obtained by filtering the original sequence c in both cases. This is the fun- 
damental property of this decomposition algorithm: even if we decompose two 
different functions, Qe^n with s G -Ei, we have to filter only one data vector to 
obtain the new set of scaling coefficients {c^ : e G Ei} and shearlet coefficients 
K,^ : £G-Bi,7er*}. 
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In the next step, the sequences and the associated functions q{e,ri),n-i are 
decomposed in precisely the same way, making use of Corollary |7.4| again. Like 
above, we filter cq twice to obtain new, further downsampled sequences C(o.o) ^^^^d 
C(o,i) together with the respective shearlet coefficients rf(o,o),7i 7 G Tq and (i(o,i),-Y, 
7 e r|. In exactly the same way we obtain C(i.o) and C(i,i) as well as (i(i,o),7, 
7 € Fq and i) ,y, 7 G by filtering ci. These first two steps of decomposition 
are illustrated in Figure [7] It can already be seen from Figure |7] that - hke the 




Figure 7. The binary tree associated with the fast shearlet decomposition. 



subdivision scheme - the shearlet decomposition becomes a binary tree labeled by 
the directional indices e. Indeed, in general we obtain the new coefficients by the 
following simple filtering. 

Algorithm 7.6. Let for some e € E be given. Then the next level of scaling 
and shearlet coefficients are computed as 

Eventually, this process ends up with coarsest level scaling coefficients c^, e G En, 
and shearlet coefficients d^^^, e G Ek, k < n, ^ Cz F*^ which describe the deviation 
from the coarse data. 

Indeed, it is now easily seen that such a decomposition must recognize "sheared" 



and thus directional components of two dimensional data since (7.2) relates, for 
e G E, the data Dw^c with the function and the respective shearlet coefficients 
must be large where the prediction by the subdivision scheme is inaccurate, i.e., at 
directional singularities. Thus, the "recipe" is to consider the shearlet coefficients 

dpkcj^ A: = 1, . . . ,n, e G £:„, 7 G F*^. 

A precise analysis of this nevertheless fundamental aspect of directional edge detec- 
tion is beyond the scope of this paper where we just want to give the framework for 
adaptive directional detections. It should also be clear that the adaptive directional 
approach is not tied to interpolatory schemes, in fact, any perfect reconstruction 
filter bank can be used as long as the projection and its complement can be ex- 
pressed properly. We plan to address these questions as well as the numerical 
implementations in a further paper, however. 
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